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1 . 

1.1 


INTRODUCTION 


The ALMA code 

The ALMA ( Axi-symmetrical Lateral Momentum Analyser) 
computer program is concerned with the numerical computation 
of two-dimensional coaxial jets with large lateral pressure 
gradients. The jets may be free or confined, laminar or 
turbulent, reacting or non-reacting. 

In its present version, ALMA predicts the mass fractions 
of the products of combustion of hydrogen and oxygen, 
based on the assumption of chemical equilibrium. 

1 . 2 Connections with previous work 

Under a previous contract (NASl-12167 ) , a computer code, 
CHARNAL , was developed for predicting the mixing and 
combustion of turbulent coaxial hydrogen-air jets. CHARNAL 
was based on the assumption that the pressure of the gases 
varied with axial distance only; lateral pressure gradients 
were regarded as negligible. ALMA accounts for these lateral 
pressure gradients. 

The common features between ALMA and CHARNAL are: 

• Both codes account for the effects of temperature and 
composition on local fluid properties such as laminar 
viscosity and specific heat of the gas mixture. 

• Both codes allow for incorporating any inlet-flow 
conditions (temperature and velocity etc.). 

• Both codes have provision for printing some or all of the 
dependent variables at pre specified axial distances, or 
after a certain number of marching steps. 

1 . 3 Layout of the report 

The remainder of this report is divided into four sections. 
Section 2 is concerned with the mathematical formulation and 
physical models employed in ALMA. Section 3 outlines the 



SIMPLE algorithm used to link the lateral momentum equation 
with both the continuity and streamwise momentum equations. 

In Section 4, sample predictions for some test cases 
specified by NASA Langley are discussed. Suggestions for 
further extensions and refinements of ALMA are included in 
Section 5. Appendix C details the input quantities required 
for using ALMA. 

2. THE MATHEMATICAL AND PHYSICAL ANALYSIS 

2 . 1 Introduction 

This chapter details the mathematical and physical basis of 
the ALMA code. Section 2.2 is concerned with the conservation 
equations of momentum, stagnation enthalpy and chemical 
species. Section 2.3 describes the turbulence model; and 
in Section 2.4 the chemical-equilibrium combustion model is 
discussed . 


2 . 2 The mean-flow conservation equations 

The differential equations listed in this section are the 
ones which express the conservation of momentum, mass, 
energy and chemical species in a two-dimensional, axi-symmetr ic, 

a|c 

boundary-layer flow with large gradients of lateral 
pressure . 


The streamwise momentum equation : 


pu 


du 

3x 
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( 1 ) 


The lateral momentum equation: (viscous terms being neglected) 
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( 2 ) 


* "Boundary-layer", in this context, refers to the 

characteristic absence of influences from downstream 
on upstream events. 
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The continuity equation : 


h. (pu) + 1 -b (pvr) = 0 


( 3 ) 


Conservation equation of th e mixture fraction ( = the 
mass fraction of hydrogen in any form) (f) : 


9f , 9f 1 9 9f. 
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Conservation of stagnation enthalpy (h) : 
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The temperature of the mixture, T, is obtained from known 
values of h,v,u and the mass fractions of the chemical 
constituents of the mixture as follows: 


A 2 , 2 
£ u +v 
h - — ^ — 


k = h = . m . h . 

J J J 


( 6 ) 


The species enthalpy ? h is defined as : 

J 


/ 


ref 


c . d T + H . m 
P,J J,T 


ref 


(7) 
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where T^ ef is the temperature of the state from which 
enthalpy is measured, is usually taken as 298.15 K, and 


H 


j,T 


ref 


denotes the heat of formation of the species 


j at T j. Equation (2.2-7) can then be written as 

T„ 


h . = 
3 


I 


c . dT + H . 000 , c „ 
p, j j ,298.15 K 


( 8 ) 


298.15 K 
T/ 


J 


0 K 


C p,j dT + ^ H j , 298 . 15 K ^ h 298.15 K h 0 


= C - T + H 
P , 3 o 


(9) 


here c , defined as 1 / c -dT , denotes a mean 

P ? J rp J P 7 J 

0 K 


specific heat, and H i s defined by: 


H o " H j , 298 . 15 K ( ' h 298.15 K h 0 K ) (10) 

2 . 3 The turbulence model 

2 . 3 .I The turbulent viscosity 

The effective turbulent-transport coefficients 

are determined by means of the k-e model of turbulence 
(ref. 1). According to this model, the magnitude of the 
viscosity depends only on the local values of the turbulence 
kinetic energy, k, on the dissipation rate of turbulence 
energy, e, and on the fluid density. This turbulent viscosity 
y t , is given by : 


y t = C y p k Ve 


( 11 ) 


where C is a constant. 
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for each point in the 
of transport equations: 


(|i) 2 - pc, (12) 
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(13) 


2.3.2 The concentration fluctuations 

Another important turbulence quantity is the time-mean square 
of the fluctuations of mixture fraction g. This quantity can 
be used to account for the effect of turbulence on the rate 
of chemical reaction (Kef. 2) or to estimate the width of a 
turbulent-flame brush (Ref .3). In ALMA, this width, which is 
defined as the region of flow over which the stoichiometric 
conditions occur, is calculated as follows. 


First, the magnitude of the mean-square fluctuations of the 
mixture fraction , g, is found from the solution of the follow- 
ing transport equation developed by Spalding (Ref. 4): 
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3x 


pv 
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JL (1 

3r ' a 


g 


is 

3r 


)+ c 


gl 


_Df\2 
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In the above equation, the 

and C n are 0.7, 2.8 and 2 
g2 


values of the constants 
respectively . 


(14) 

o , C i 
g gl 


Secondly, a battlement-shaped distribution of the mixture 
fraction vs time is assumed (Ref. 3); this means that the 
instantaneous mixture fraction fluctuates rapidly between 
an upper and lower value, and spends no time in between. 
From this distribution, and from the values of f and g, the 
upper and lower values, f , and f , of the instantaneous 

~r — 

mixture fraction are calculated. 
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Finally, the width of the turbulent flame brush at any 

downstream location is obtained by determining the two. 

radii at which f , and f equal f 

+ - st 


2 . 4 The combustion model 
2,4. l Equat ions 

Reference (5) describes an equilibrium-chemistry model, 
developed and used in conjunction with the CHARNAL computer 
program to predict the properties in a hydrogen-oxygen 
flame. The main features of the model are described below. 

Four equilibrium reactions are assumed as follows: 


0 + 0 -+ 0 2 (15) 

H + H -+ H 2 (16) 

H + OH -+ H 2 0 (17) 

0 + H -> OH (18) 


The six species involved in these reactions are considered 
to be present in the mixture, together with nitrogen, which 
is inert. To assist in developing the equations to predict 
the equilibrium concentration of the species, two intermediate 
quantities are defined, namely: 


W 




W 


f = m tI + m TT + 


m. 


H 


m. 


'II 2 • “H "H 2 o w oh oh 


w w 

= m 0 o + m o + WITT m OH + w m 

Z UH 


H 2 0 H 2 0 


(19) 


and F 


+ 


+ 


( 20 ) 



where f is the total mass fraction of- hydrogen in any form, 
and F is the total mass fraction of oxygen in any form. 

Also required is the auxiliary relation which states that 
the sum of the mass fractions of all elemental species must 
be unity: 

f + F + m„ = 1 (21) 

One assumption is now introduced: it is presumed that the 

coefficients of turbulent diffusion of the chemical species 
are equal to each other at every point in the flow. A well- 
known consequence is that f is linearly related to m w ; and 
the constants in the relation can be determined from the 
boundary conditions: 


m 


m 


N 


N 


2 

2 , air 


1 


f 


( 22 ) 


What remains now is to define the necessary equilibrium 
constants . 


As 



is common in thermodynamics, the equilibrium 
is defined for the reaction aA + bB > cC 


constant , 
by : 



X~ c-a-b 

_C P 

y a v b 

X A X B 


(23) 


where p is the pressure in atmospheres and X . is the mole 

J 

fraction of species j. For each of the four reactions in 
the present model, (c-a-b) equals -1. 

It is convenient to express the concentrations in terms of 
mass fractions using the relation: 


I 


7 


( 24 ) 


m . 
J 


W . 
= _A 
W 


X. 

J 


Substitution of equation (24) into (23) permits derivation 
of the modified equilibrium constant K^, thus: 



W 

c 

K p W — t- 

P \y a 

"B 


m. 


a b 
m A m B 


(25) 


The modified equilibrium constants for the reactions (15), 
(16), (17) and (18) can then be written as: 
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H, 


k; = 
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0 f 


k: = 
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m H 0 
K 3 “ ETra 


H OH 


and 



m 


OH 


m 0 m H 


(26) 


(27) 


(28) 


(29) 


If thermodynamic equilibrium prevails, the K"'s take values 
which depend upon temperature alone. Their values are given 
in Appendix B. 

In the above set of equations, there are 7 unknowns 
(m R , m Q , m Q , m QH , m H Q , m N and F; f being given) and 
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7 equations { (19), (20), (21), (26), (27), (28) and (29)} . 

The problem is therefore soluble. 


2.4.2 Solution procedure 

The solution procedure is an iterative one, operating on the 
two equations: 



(31) 


Here the quantities A.B,C,D are defined by: 

A - — (32) 





(35) 
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These quantities depend upon temperature alone, and can be 
tabulated at the start. 

Equations (30) and (31) have been derived from equations 
(19), (20), (21), (26), (27), (28) and (29) by straight- 
forward algebraic manipulation, including the use of the 
formula for the roots of a quadratic equation. 

Because the RHS of equations (30) and (31) contains and 

m ¥T as well as f and F, iteration is necessary. Direct 
H 2 

substitution of the latest guesses for m„ and m^ is 
employed in ALMA, and has been found to be satisfactory in 
the circumstances which have been investigated; but a 
Newton-Raphson technique could be employed if divergence 
were encountered. 

3. THE NUMERICAL SOLUTION PROCEDURE 
3 . 1 Introduction 

The numerical procedure, employed here to solve the set of 
partial differential equations that govern the flow and 
combustion of axi-symmetr ic boundary layers, is of the 
marching-integration type. Because of the need in supersonic 
flows to solve the lateral momentum equation (2) 
simultaneously with the other conservation equations, a new 
algorithm (SIMPLE = Semi implicit Pressure Linked Equations) 
is introduced to link the streamwise and radial momentum 
equations with the continuity equation, so as to obtain a 
pressure field. The remainder of this chapter outlines the 
basic features of this algorithm and of the numerical 
procedure . 
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3.2 The finite-difference grid 

shows a portion of the f inite— difference grid on 

the x-r and the x-w coordinates, where w is a dimensionless 
stream function (Ref. 6). 

The variables u(the x-wise velocity component), p(the 
pressure), and p(the denisty) are all defined at the grid 
nodes. The lateral velocity component, v, is defined at 
locations (o' s) midway between those of the u-velocity, as 
shown in Figure 1, in conformity with usual "staggered-gr id" 
practices. The definition entails that the pressures, of 
which the differences affect v, straddle the v location. 

3 . 3 The formulation of the solution algorithm 

Three main steps are followed here to formulate the solution 
algorithm. T he first step is to use the Patankar-Spald ing 
(Ref. 6) procedure to obtain a formula which connects the 
value of a dependent variable at a downstream node with its 
value at an upstream node. The second is to use the lateral 
momentum equation to derive a formula which relates the 
radial velocity component at a node (i + £) to the pressures 
at the nodes (i + 1 and i) which straddle the v-node at the 
same axial location. The third step is to derive a"pressure 
correction" equation and use it to determine a pressure field 
which satisfies both the streamwise and lateral momentum 
equations and the continuity equation. 

Each of the above three steps is discussed below. 
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3.3.1 The Patankar-Spalding procedure 

This procedure is detailed in Reference 6, therefore only the 
main features are discussed here. 

(i) The primary differential equations ( 1 ) , ( 4 ) , ( 5) , ( 12 ) , 
(13) and (14) are transformed so that the independent 
variables are the longitudinal distance x and the 
dimensionless stream function co defined as: 

m = O' - f j ) / ( g - Vj) (36) 

where ¥ is obtained from 


ay 

Sr 


pur 


Subscripts I and E refer respectively to the "inner" 
and "outer" boundaries of the flow. The resulting 
differential equations are all of the form: 


3x 


+ (a + bo3 ) 


30 _ 3 
3o) 3co 


-i> + d 


(37) 


where the terms on the left-hand side represent 
convection and those on the right hand side express 
the diffusion and source of the variable 0. The 
coefficients a and b are functions of the entrainment 
rates, whilst c involves the effective diffusion 
coefficient . 


12 



(ii) The differential equations are expressed as finite- 
difference equations, connecting the values of 0 at 
the nodes i, i + 1, i - 1 (Figure 1) at a downstream 
station D with those at the corresponding nodes at the 
upstream station U. The general form of the finite- 
difference equation is: 


D. <A . _ = A . d> . . ^ _ + B. <& . ^ _ + E . d> • TT + G . d> . , - TT 
l *i,D l *1+1, D l y i-l ,D l *i,U l *1+1, U 


+ H 


• $ • 1 tt + S . 

1 * 1 - 1,11 1 


(38) 


where A., B., etc., are treated as constants, the 

l l ’ 

expressions for which are deduced by integration of the 
differential equation (37) over a control volume 
surrounding the node where prevails. 


(iii) The integration proceeds by "marching” downstream, the 
values of (jt> at station D being calculated from those 
at station U. 


At each forward step, new values are ascribed to 
and H'j, the stream functions at the grid boundaries. 
These values, together with the continuity equation, 
determine the geometrical locations of the boundaries. 


3.3.2 The relation between v and p 

The relation between v and p is obtained by integrating the 
lateral momentum equation (2). The finite-difference form 
of this relation is: 


V i+*,D V i+£,U 


5 x 


" p i- p i+l 

pu 

i+i ,D 

r i+l~ r i 


(39) 


where the notations of Figure 1 are used. 
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Equation (39) shows that a solution for v^ + ^ D cannot be 
obtained in an explicit manner since the radial pressure 
distribution at D is not a pre-specif ied quantity. This 
problem is dealt with in the following section. 

3.3.3 The pressure correction procedure (SIMPLE) 

The central idea of the pressure-correction procedure is to 
make an informed guess of the pressure distribution at the 
downstream station D; then corrections are added to the 
pressures, and to the velocity components u, v, to ensure 
compatibility of flow areas. 


In this section it is aimed 
equation of the form: 


" C i + i P i + 1 


- c . i p . . 

i i-l 


to derive a pressure-correction 


+ 



+ 


c 


i-i 






(40) 


where the superscripts (*) and {') denote respectively a 

guessed quantity and a correction component. In equation 

(40) the right-hand side is the area-discrepancy appropriate 

* 

to the guessed velocity field; A being the flow area 

associated with ,the continuity and streamwise momentum 
* 

equations and a being the area associated with the lateral 
momentum equation. The c's are coefficients having the 
dimensions of area divided by pressure. 


It is important to note that the solution of the pressure- 
correction equation (40) yields a pressure-correction 
component which, when added to the guessed pressure, cancels 
the discrepancy in the flow areas. 


Equation (40) will now be derived. 
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Three steps are followed in the derivation. The first is to 
obtain a p^a"* relation; the second is to find ap^A' 
relation, and the third is to combine those two relations to 
produce equation (40). 

p^a.' relation 

It is sought to obtain an equation in the form: 


a 


= c 


i+i 


(j»I - 


p LP 


+ c 


i-4 


( p± - 


p i-i> 


(41) 


where c's are coefficients having the dimensions of area 
divided by pressure. The starting equation is that which 
links v and p, i.e. equation (39). 

Let the downstream quantities u , v , p , p , A, a , r be composed of 
two components according to: 


p = 

p* 

+ 

P 

(42) 

u = 

u* 

+ 

u" 

(43) 

v -- 

v* 

+ 

* 

V 

(44) 

p 

p* 

+ 

p 

(45) 

A = 

A* 

+ 

A' 

(46) 

a = 

a* 

+ 

a.' 

(47) 

r = 

r* 

+ 

x' 

(48) 


Substitution into equation (39) of the values of p,u,v,p,r 
from equations (42), (43), (44), (45) and (48), followed by 
subtraction of the resulting equation from (39) yields: 


v 


i+4 


(J* ) 

v p*u* ’ i+i 



(49) 


where p, u and r are approximated by their guessed values. 



Another equation is required to connect a 7 with v / ; this 
is obtained by referring to Figure 1. The cross-stream 
area a . , at D, for the shaded control volume associated with 
node i can be calculated from: 


a . 

1 


a . TT + 

l ,U 


J r 

L i+l,U 


tan a 


i + 4 


r . i TT tan a . , 1 6 x 

i-4, u i-4 r 


(50) 


where a. T , is the upstream cross-stream area. The 
i,U 

inclination a of the constant-w lines is given by: 


a 


• 1 1 

v _ m_ 
u pu 


(51) 


• H 

where m is the rate of mass flow across the constant-a) line. 


Substitution into equation (50) of the values of a,r,v,u,v' 
from equations ( 51 ) , ( 48 ) , ( 44 ) , ( 43 ) and (49), followed by 
subtraction of the resulting equation from equation (50), 

leads to the desired equation (41), where c now has the 

2 2 
value r 5 x/(6r p* u* ). 

p^A^relation 

The desired p"vA''equat ion is 



which expresses the effect of pressure change on the down- 
stream cross-stream area. Three equations are employed to 
derive equation (52). These are the continuity, the 
streamwise momentum equations and the gas-law relation for 
isentropic expansion. 

The mass flow rate across the downstream face of the shaded 
control volume (Figure 1) is: 



( 53 ) 


p. u. A. = (f - ! )(d) . - a) . ,) 
1 l l v E I' v i+i l-i 


* * * 
= u. A, 


Substitution into equation (53) of the values of p,u,A from 


equations ( 45 ) , 

(43), 

(46) gives 

* 
P • 

u: 

a: 

i 

+ — | 

+ ~ = 0 

T’ 

P i 

u. 

1 

A. 

l 


where the terms containing the cross products p^u^ p^A^ u^A^ 
p^u^A^ have been neglected. 


Now the variation of u with pressure can be estimated by 
noting that, if the shear stresses and body forces are 
supposed to be uninfluenced by the pressure change, the 
streamwise momentum equation (1) leads to: 


du 1 

dp pu 


(55) 


The substitution of equations (42), 
yields : 


u . 

l 



(43) , (45) into (55) 


(56) 


Further if the process of pressure 
be isentropic, then: 



adjustment is assumed to 


(57) 


Now the p' ^ Adequation (52) can be easily deduced by 
combining equations (54), (56) and (57). 
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Finally the pressure-correction equation (40) can now be 

obtained by noting that the net effect of the pressure- 

correction process is to cancel the discrepancy in the 

calculated cross-stream flow areas, i.e. A^ should be equal 

to a. or 
l 


A. + A'. = a. + a' 

1111 


(58) 


The direct substitution of equations (41), (52) into (58) 

produces the pressure-correction equation (40) for the node 
i; similar equations are obtained for all the nodes across 
the grid. 


3 . 4 Conditions at the E-boundary in superson ic f l ow 

When the free stream is supersonic, the values of u and v 
at the E-boundary are evaluated from the following two 
equations : 


and 


v 2 V 1 


U 2 U 1 


P 2 -Pf 


(P 2 ~ P 1 )/(P 1 u x ) 


Pl U l 


tan 


[a +sin 1 ( )1 - 

L 1 l/~2 2 

/ u l +v i 


tan a. 


(59) 


(60) 


where the notations of Figure 2 are used, and c 1 is the speed 
of sound at station 1. 
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Equation (59) is derived from the one-dimensional streamwise 
momentum equation for the shaded control volume (Figure 2). 
Equation (60) is obtained by combining the equations of 
continuity and lateral momentum for the shaded control 
volume and then substituting for the angle (0 -a^) between 
the streamline and the shock from small-wave theory. f 

It is seen from equation (60) that the velocity v^ at the 
E-boundary can be calculated in terms of quantities 
immediately upstream. The pressure a t the E-boundary is 
set equal to the value calculated at the near-boundary node. 

3 . 5 The basic solution steps 

The solution procedure will now be reviewed, based on a 
consideration of the portion of the finite-difference grid 
shown in Figure 1. The variables at station D are to be 
determined by the integration of the governing differential 
equations over the shaded control volume. This integration 
is performed as follows: 

a) The pressure gradient across the forward step is guessed; 

* 

the pressure p at station D is then calculated. 

b) Equation (1) is integrated as in Ref. (6). The resulting 
velocities at D are only estimates since they result 
from a guessed pressure field. 


\ From small-wave theory, the angle (9 -a^) is given by 

1 


6 - a. = sin 1 [rj-j] 


K +V l 


where c^ is the speed of sound at station 1 
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c) The continuity equation (3) is integrated to yield the 

radial locations (r.) of the constant w lines at D. The 

1 

area A. of the downstream face of the control volume 
1 

(shown shaded in Figure 1) is deduced therefrom. 

d) The lateral momentum equation (2) is integrated to 

* 

yield the v velocities at D; these velocities are only 
estimates, since the pressure field is a guessed one. 

e) The inclinations (a^) of the constant-ai lines are 

calculated from the values of the v velocities at D. 

Now it is possible to re-evaluate the area of the 

downstream face of the control volume from the knowledge 

of a.. The new area value is denoted a., 
l i 

f) The values of A. and a. would be identical if the 

ii 

guessed values of the pressure field were the correct 
ones. In practice, however, there is a discrepancy 
between A. and a. : for the former was calculated from 

l l 

the streamwise momentum and the continuity equations, 
whereas the latter was calculated from the lateral 
momentum equation. 

g) In order to eliminate the discrepancy between the values 

of a^ and A^, a pressure-correction component is 
introduced into the analysis with the net effect of 
yielding a correct pressure field that will satisfy 
simultaneously the equations of: streamwise momentum, 

lateral momentum and continuity and hence cancels the 
imbalance in the values of a. and A.. The pressure 

ii 

corrections are obtained from solving the set of equations 
such as (40) by the tri-diagonal matrix algorithm. 

h) The guessed pressure field and the associated velocity 
fields are then augmented by the calculated correction 
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components; the result is a more correct flow field, 
which certainly satisfies continuity, and which satisfies 
momentum as accurately as is desired . 

i) The above operations (a to h) complete one integration 
step; they are to be repeated for successive forward 
steps until the whole domain of integration has been 
covered. 

4. DISCUSSION OF SAMPLE PREDICTIONS 

To demonstrate the general capabilities of ALMA, five test 
cases have been computed; these involve the mixing and 
combustion of two coaxial jets of hydrogen and air. The 
geometry and the initial conditions for the test cases are 
given in Figure 3 and Table 1 . 

In this section, the results are presented in the form of 
radial profiles of the quantities u , v , T , k , e , f , g , £ and of the 
mass fractions of the species 0^, E^, H^O, H, 0, OH and Ng 
in Figures 4 to 17 inclusive. Because the values of the 
above quantities vary by orders of magnitude and are plotted, 
at a given station, for convenience in two figures, each is 
normalized by reference to its maximum and minimum values in 
the profile, so that the normalized values extend from 0 to 1 . 

4 . 1 The j e t surrounding a centre-body 

Case 1 . The results are discussed by examining the radial 

distributions of the quantities u , v , T , k , e , f , g and £ and the 

mass fractions of species at three different downstream 

stations: x/Y . = 20, 60 and 100 (Figures 4, 6, 8, 5, 7 and 9). 

3 

Profiles at x/Y^ = 20 

Figure 4 shows that at x/Y. = 20 (i.e. nearest to the nozzle) 

3 

all the predicted quantities behave as expected. For 
example, the u-profile exhibits a steep gradient near the 
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Table 1. Initial Conditions of the Test Cases 


Test 

Case 

1 

2 

3 

4 

5 

p i 

C 10 5 N /nvj 

1. 

1. 

1. 

2. 

1. 


(10 5 N/m 2 ; 

1 . 

1. 

1. 

1. 

1. 

T i 

( K) 

316 . 

180. 

247. 

247. 

247. 

t e 

( K) 

885. 

885. 

1459. 

1459 . 

1459. 

T w 

( K) 

300. 

300. 

- 

- 

900. 

U I 

(m/s ) 

1303. 

2100. 

2389. 

2389. 

2389. 

U E 

(m/s) 

1800. 

1800. 

1432. 

1432 . 

1432. 

Y . 

(m) 

1. 00264 

0 . 00264 




J 






R . 
J 

(m) 

- 

- 

0. 00238 

0.00238 

0.00238 

R 

(m) 

j. 01214 

0. 01214 

— 

— 

0. 01633 

w 






in 

°2 

) outer 
) 

.232 

. 232 

. 232 

. 232 

.232 

m 

N 2 

) jet 
) 

. 768 

. 768 

. 76 8 

. 76 8 

. 76 8 

m 

«2 

) inner 
) jet 

1. 

1. 

1. 

1. 

1. 
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centre-body wall and a less steep one in the shear layer 
between the jet and the free stream. 

Also, as is well known, the k-profile (Figure 4) attains a 
maximum at about the middle of the shear layer where both 
u and v have large gradients. The dissipation rate of 
turbulence energy is highest near the wall, where the length 
scale of turbulence is at its minimum. 

Chemical reaction is evident by the peak in the temperature 
profile ^Figure 4) at the radial location where the mixture 
fraction, f, attains its stochiometr ic value, f g ^_, of about 
0.028. 

The mean square of f-f luctuat ions , g, reaches its maximum 
value where f has its steepest radial gradient; this is as 
expected . 

Figure 5 shows that the radicles 0 and OH have their maximum 
concentrations at the peak temperature, whereas H attains 
its highest value at a slightly lower temperature towards 
the hydrogen-rich region. 

It is also displayed in Figure 5 that the profiles of and 
Og overlap in a small region centered around the location of 
f gt . The two species N 2 and H 2 0 exist in considerable 
concentration over the whole domain, thought their peaks lie 
in the air-rich and the hydrogen-rich regions respectively. 

Profiles at x/Y . = 60 
J 

Figure 6 shows thaf at a distance 60 Y. downstream, the shear 

J 

flow develops and the initial step-shaped u-profile becomes 
relatively smooth. The maximum values of k and e are 
consequently lower than those at the upstream station of 20 Y \ . 
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The temperature profile still has its maximum near the 

location of f ^ - The value of this maximum is 30% higher 

than that at the station of 20 Y.; this is due to the 

3 

development of chemical-reaction as indicated by the lower 
value of iiUj (Figure 7). 


Profiles at x/Y. = 100 
' J 


Figure 8 shows that at a distance of 100 Y., the u-profile 

J 

is similar to that at 60 Y ^ . Also, the profiles of k and 
£ are similar to those at 60 Y^ , though the level of k is 
slightly lower and the maximum £ is higher than that at 60 


Y . . 
J 


Figure 9 displays the profiles of the species concentrations 

at x/Y. = 100. Here the maximum value of m^ , near the wall 
J Hg 

•of the centre-body, is only 15.8% of its initial value at 
the nozzle exit. 


Comparison between Cases 1 and 2 . Now the predictions of 

case 2 at 100 Y. (Figures 10 and 11) are compared with those 
3 

of case 1 at the same station. Both cases have the same 
geometry and initial conditions, except that in case 2 the 
initial u-velocity of the hydrogen jet is 50% larger than 
that in case 1. 


Figure 10 shows that the u-profile is not yet fully developed; 

it still retains the step shape of the initial distribution. 

This behaviour is explained in what follows. Here, the 

temperature level is lower than in case 1, because the shear 

layer is predominatly richer in fuel than in case 1; the 

maximum values of f and m u are 0.31 and 0.286 as compared 

z 

to their respective values of 0.182 and 0.158 in case 1. 

Lower temperatures result in higher densities and hence a 
lower rate of velocity-decay. 



4 . 2 The jet without a centre-body 


Case 3 . Figures 12 and 14 show the radial profiles of u,T, 

k,f,g and £ at two downstream-stations: x/D. = 10 and 50. 

0 

The corresponding profiles of the species mass-fractions are 
shown in Figures 13 and 15. 

At a distance of 10 downstream (Figure 12) , k reaches its 
maximum in the high-shear region away from the axis, whereas 
at a distance of 50 (Figure 14) where the shear flow is 
developed, k attains its peak at the axis. This behaviour 
is generally a feature of axi-symmetr ic jets and is different 
from that of the jet in cases 1 and 2 because of the presence 
of the centre-body. 

The temperature profiles at both axial locations have their 
highest value where f equals f ; this is also where the 
radicles O and OH reach their maximum concentrations 
(Figures 13 and 15). 

Comp arison between Cases 3 and 4 . The profiles of u,T,k,f, 
g and £ at a downstream distance of x/D^ = 50 for case 4, are 
shown in Figure 16. 

The profiles of the species mass-fractions at the same 
distance are shown in Figure 17. 

The initial conditions of case 4 differ from those of case 3 
only in that the static pressure of the hydrogren jet of the 
former is twice that of the latter. The effect of this 
higher pressure is seen by comparing Figures 14 and 16. It 
is evident that both u and f have a lower rate of decay in 
case 4 than in case 3; this can be attributed to the higher 
mass flow rate in case 4 caused by the higher pressure and 
density. Also the maximum mass fraction of at a distance 
of 50 D. in case 4 (Figure 17) is higher than that of case 3. 

J 
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4 . 3 The confined jet, Case 5 . 

Regarding case 5, the confined-flow case, it was planned to 
continue the integration to a downstream position of 50 D.: 
however, it was not possible to proceed with the computation 
farther than 10 , because choking was predicted there. 

This choking is detected by the computer program when the 
calculated flow-area necessary to satisfy the conservation 
of lateral momentum exceeds that required for conservation 
of mass (here this is equal to the area of the confining 
duct). This situation indicates that the intitial flow 
conditions for case 5 are not compatible with the flow 
through the duct. 


It is suggested that new initial conditions should be 
specified for this case, e.g. by successively lowering 
the velocity of the hydrogen jet, and performing further 
runs . 
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Figur e 1: The finite-difference grid in both the x-q 3 

and the x-r coordinates. 




r - Shock 





c - Case 5 


Figure 3 


: The Flow Geometry of the Five Test Cases 
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Fig. 8 Radial profiles of u , T , k , f , £ , v at x/Y_.=100, Case 1. 
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Radial profiles of the species mass-fractions at 
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Fig. 17 Radial profiles of the species mass-fractions at 
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5 * SUGG E STION S FOR FURTHER EXTENSIONS AND REFINEMENTS OF ALMA 
The following extensions and refinements to ALMA are easily 
possible; and they may prove to be necessary when the 
predictions of ALMA are compared with experimental data:- 

C a ) I nf ro ducti on of realistic chemistry 

As has been explained, thermodynamic equilibrium has 
been postulated in ALMA; the chemical reaction rates 
have been presumed to be very fast. The presumption 
is at best an approximation to the truth; and it 
will undoubtedly be necessary, from time to time, to 
examine how the predictions will change if chemical- 
kinetic limitations are allowed for. 

These limitations may be of two kinds; those associated 
with the chemical kinetics of laminar-flow processes; 
or those concerned with "eddy break up". New ideas 
(Ref. 2) concerning the interactions of these two modes 
are now available; and the ALMA code can be easily 
extended to accommodate these ideas. 

The extension naturally takes the form of the provision 
and solution of further differential equations, together 
with auxiliary relations concerning sources. However, 
in order that the associated computer time does not 
become excessive, simultaneous improvements in the 
numerical procedures discovered since ALMA was conceived, 
may be advantageously introduced. 

( b ) Introduction of "part i al ly-parabolic" effects 

ALMA is an advance over CHARNAL, in that it takes account 
of the lateral-momentum effects in the supersonic-flow 
regions. However, where the flow is subsonic, the variation 
of pressure with radius has to be ignored if the marching 
nature of the integration is to be preserved. 
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In some circumstances , the inaccuracies resulting 
from this practice are not to be tolerated. 

Fortunately, experience now exists of handling the 
lateral-momentum effects in subsonic regions by the 
so-called "part ially-parabolic" method (Ref. 7 ). 

It is true that repeated marching integrations must be 
made, to this end; so computer time increases. 

However, the partially-parabolic method is economical 
of storage: only the pressure has to be stored for the 

whole two-dimensional flow field; so the storage 
requirement is scarcely greater than that of ALMA. 

( c) Allowance for small patches of recirculation . 

In a still further extension, "elliptic" flow regions 
may be incorporated. These require still further 
storage. However, recent coding developments (Ref. 8' ) 
have permitted the additional storage, and the more 
numerous iterations, to be restricted to the regions of 
recirculation themselves. It is not necessary, just 
because shock-wave-boundary- layer separation occurs in 
one area, to forego entirely the economies associated 
with marching integration. 

This refinement will also be important if it should 
become important to compute the influence of the finite 
thickness of the "lip" of the hydrogen-injection nozzle. 

Of course, the stream-line coordinate system must be 
abandoned, at least in the recirculation region. 

However, means exist of retaining its advantages 
elsewhere. 
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APPENDIX A. Symbols 


C ,C, ,C 0 ,C ,C _ 

\i 1 2’ gl g2 


P,3 


D . 
J 

F 

f 


st 


g 

h 


h j,t 


ref 


K 

P 

K. 

i 

k: 

i 

k 

l 




cross-sectional flow area of a control 
volume in the continuity equation; 

cross-sectional flow area of a control 
volume in the lateral momentum equation 

coefficients in the pressure correction 
equation ; 

constants coefficients appearing in 
turbulence model; 

speed of sound; 

mean specific heat of species j; 

diameter of jet nozzle; 

mixture fraction of elemental oxygen; 

mixture fraction or mass fraction of 
elemental hydrogen; 

stoichiometric value of f; 

mean square fluctuations of f (= f' 2 ): 

enthalpy ; 

heat of formation of species j at T = T 

stagnation enthalpy; 
general equilibrium constant; 
equilibrium constant for the reaction i 
modified equilibrium constant ; 
kinetic energy of turbulence; 
length scale of turbulence = k 3//2 /e; 
mass fraction of species j 


ref 
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p 


R. 

J 


R 


w 


T 

u 


V 

W 

W . 
3 

x 

X . 
J 

y 


Y . 
3 


static pressure; 

radius (measured from axis of jet); 

radius of hydrogen pipe at exit; 

radius of duct enclosing the jets or 
radius of center-body; 

Temperature ; 

axial component of velocity (along a 
constant-w line); 

radial component of velocity; 

average molecular weight of the mixture 

molecular weight of j species; 

distance along axis of jet; 

mole fraction of species j ; 

radial distance (measured from axis of 
jet); 

radius of the outer edge of the jet 

width of the annulus of the hydrogen 
jet ; 


46 



Greek Symbols 


a 


Y 

0 


e 

P 

0 

Y 

0 ) 

Superscripts 


Symbol 

angle between the constant-u line and 
axis of jet; 

the ratio of specific heats; 

angle between a streamline and axis of 
jet ; 

rate of turbulence energy dissipation; 
density; 

general variable; 
stream function; 

non-dimensional stream function; 


Symbol 

estimated quantity; 

correction component except in the 
definitions of the chemical equilibrium 
constants . 
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APPENDIX B. Table of Equilibrium Constants 


Temp. 

k: 

i 

■5 


K 4 

200 

5.0+100 

3.1E+99 

2.6+100 

2.6+100 

400 

2.8E+52 

2.7E+58 

3.2E+51 

6.2E+59 

600 

2.3E+33 

4.3E+36 

5.8E+32 

9.0E+37 

800 

5.9E+23 

4.9E+25 

2.2E+23 

9.3E+26 

1000 

9. 7E+17 

1.2E+19 

4.5E+17 

2.2E+20 

1200 

1.2E+14 

5.0E+14 

7.0E+13 

8.2E+15 

1400 

2.1E+11 

3.5E+11 

1.2E+11 

5.4E+12 

1600 

1.7E+09 

1.1E+09 

1.1E+09 

2.2E+10 

1800 

3.9E+07 

2.1E+07 

2 . 7E+07 

3.0E+08 

2000 

1.9E+06 

7.0E+05 

1.4E+06 

9.8E+06 

2200 

1.5E+05 

4. 3E+04 

1.2E+05 

5.8E+05 

2400 

1.9E+04 

4.2E+03 

1.6E+04 

5.5E+04 

2600 

3.4E+03 

5.8E+02 

2.8E+03 

7.6E+03 

2800 

7.5E+02 

1-0E+02 

6.4E+02 

1.3E+03 

3000 

2.0E+02 

2.4E+01 

1.74+02 

3.1E+02 

3200 

6.3E+01 

6.8E+00 

5.7E+01 

8.5E+01 

3400 

2.3E+01 

2. IE +00 

2.1E+01 

2.7E+01 

3600 

9.3E+01 

7.9E-01 

8.6E+00 

9.8E+00 

3800 

4.1E+00 

3.2E-01 

3.8E+00 

3.9E+00 

4000 

1.9E+00 

1.4E-01 

1.8E+00 

1.7E+00 

4200 

1.0E+00 

6.8E-02 

9.8E-01 

8.3E-01 

4400 

5.5E-01 

3.5E-02 

5.4E-01 

4.2E-01 

4600 

3.2E-01 

1.8E-02 

3.1E-01 

2.2E-01 

4800 

1.9E-01 

1.0E-02 

1.9E-01 

1.3E-01 

5000 

1.2E-01 

6.4E-03 

1.2E-01 

7.7E-02 

5200 

7.9E-02 

4 o0E-03 

7.8E-02 

4.8E-02 

5400 

5.3E-02 

2o5E-03 

5. 3E-02 

3 . ] E-02 

5600 

3.6E-02 

1. 7E-03 

3.6E-02 

2.0E-02 

5800 

2.6E-02 

1.1E-03 

2 „ 6E-02 

1.4E-02 

6000 

1.8E-02 

8.1E-04 

1.8E-02 

9.8E-03 


A 

B 

C 

5 

4.4E-51 

1.7E-50 

2.42E+50 

1 . 95E+00 

5.9E-27 

6.0E-30 

4.3E+29 

2.19E+04 

2.0E-17 

4.7E-19 

1.09E+19 

8.46E+02 

1.2E-12 

1.4E-13 

4.95E+13 

1.56E+02 

1.0E-09 

2.7E-10 

2.94E+10 

5.94E+01 

8.7E-08 

4.4E-08 

2.14E+08 

3. 14E+01 

2.1E-06 

1.6E-06 

5.22E+06 

1.81E+01 

2.4E-05 

2.5E-05 

3.67E+05 

1.32E+01 

1.5E-04 

2.1E-04 

4.53E+04 

9.45E+00 

7.2E-04 

1.1E-03 

8.63E+03 

7.76E+00 

2.5E-03 

4.8E-03 

2.24E+03 

6.96E+00 

7.0E-03 

1.5E-02 

0.71E+03 

5.77E+00 

1.7E-02 

4-1E-02 

2.6E+02 

5.30E+00 

3 „ 6E-02 

9.6E-02 

1.11E+02 

4.49E+00 

7.0E-02 

2.0E-01 

5.37E+01 

4. 37E+00 

1.2E-01 

3.8E-01 

2.95E+01 

3.87E+00 

2.0E-01 

6.7E-01 

1.70E+01 

3.62E+00 

3.2E-01 

1.1E+00 

1.02E+01 

3.45E+00 

4.9E-01 

1. 7E+00 

6.38E+00 

3.25E+00 

7.0E-01 

2.6E+00 

4. 30E+00 

3.09E+00 

9.8E-01 

3.8E+00 

3.12E+00 

3.09E+00 

1.3E+00 

5.3E+00 

2.20E+00 

2.89E+00 

1.7E+00 

7.2E+00 

1.59E+00 

2.69E+00 

2.2E+00 

9.6E+00 

1.30E+00 

2.79E+00 

2 . 8E+00 

1.2E+01 

9.62E-01 

2.59E+00 

3 „ 5E+01 

1 o 5E+01 

7.49E-01 

2.52E+00 

4 o 3E+00 

1.9E+01 

6.20E-01 

2.53E+00 

5.2E+00 

2.4E+01 

4.85E-01 

2.49E+00 

6 . 1E+00 

2.9E+01 

4.22E-01 

2.47E+00 

7.2E+00 

3.5E+01 

3.44E-01 

2.47E+00 
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AP PENDIX C . Input Quantities Required for ALMA 


ALMA solves 7 simultaneous partial differential equations 

Oi 

for the dependent variables u,v,h,k,e,f and g and employs 
typically 21 nodes in the radial direction. The program is 
written in basic FORTRAN IV language and is ready to run 
on the CDC 6600-series computers (although only minor 
modifications are required to adapt it to the UNIVAC 1106/ 
lllOor IBM 360 machines). ALMA requires approximately 
5 seconds compilation time on a CDC 6600 (FTN compiler) and 
about 40000 decimal storage locations. 

Castellated profiles are assumed for the variation of the 
dependent variables between the grid nodes and a typical 
forward step size of 0.2 times the local shear layer width 
is employed. Approximately 5 axial steps can be executed 
per second. 

The remainder of this appendix is concerned with the 
information that the user has to provide via data cards and 
data statements in subroutines MAIN and THERM. 

Subrout i ne MAIN . 

( i ) Control indices in chapters 0 and 1 : 

Fortran symbol Significance 

ITDIM: Max. number of nodes in the radial 

direction included in the profile 
plot ( >N+1 ) . 
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UI : 

Axial velocity of the inner fluid 
(m/s) 

TI : 

Temperature of the inner fluid ( K) 

PI: 

2 

Pressure of the inner fluid (N/m ) 

02E : 

Mass fraction of oxygen in the 
outer fluid 

H20E : 

Mass fraction of water-vapor in the 
outer fluid 

AN2E : 

Mass fraction of nitrogen in the 
outer fluid 

TWALL: 

Wall, temperature ( K) 

Subroutine THERM. 


Control indices in 

entry BEGIN 

NR: 

Number of reactions 

NS: 

Number of species 

HO: 

The quantity {H j>298 15 K ~ 
<h 298.15 K -h o K)} in equ. (2.2-9) 


WM: 


Molecular weight of each species 
considered in the reaction. 



JTDIM: 

ILDIM: 

JLDIM: 
KASE: 
NDATA : 

NS TAT: 

NPROF: 

DXPROF : 

NPLOT: 
DXPLOT: 
LASTEP : 

ISTABL: 


Max. number of profiles plotted at 
any axial location. 

Max. number of stations included in 
the longitudinal plots (^LASTEP) 

Max. number of longitudinal plots. 

Test case number. 

Controls whether the input data 
of the equilibrium constants and 
enthalpy will be printed (NDATA =1) 
or not (NDATA f 1) 

Number of steps between output 
stations. 

Number of steps between profile 
stations 

Axial distance between profile 
stations 

Number of steps between plot stations 

Axial distance between plot stations 

Number of steps at which calculation 
terminates. 

Adjusts, if set equal to 2, the 
forward step DX according to local 
mach number. 
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NREACT: 


Controls whether the flow is 
chemically reacting (NREACT=2) or 
non reacting (NREACT=1) 

MODEL: Equals 2 for turbulent flow and 

equals 1 for laminar flow. 


(ii) Control indices and input data in chapter 2 . 


N: 


Number of grid intervals in the 
radial direction (N<20) 


FRA: 


Specifies forward step length as a 
fraction of the width of the shear 
layer. 


KRAD: 


Specified whether the flow is 
axi-symmetric (KRAD=2) or plane 
(KRAD=1) . 


YOUT: 


Width of jet-annulus or diameter or 
jet-nozzle (metres) 


XULAST : 


Maximum streamwise distance (metres) 


KIN: 


KEX: 


Specifies 

inner boundary condition. 

KIN * 1; 

wall boundary 

KIN = 2; 

free boundary 

KIN = 3; 

symmetry boundary 

Specifies 

outer boundary condition 

KEX = 1; 

wall boundary 

KEX = 2; 

free boundary 

KEX = 3; 

symmetry boundary 
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KONFIN: Specifies flow confinement 

KONFIN = 1; confined flow 
KONFIN = 2; free flow 

CSALFA: Cosine of angle of inclination of 

streamlines to axis of symmetry 

DXLIM: Controls forward step size so that 

entrainment per step is limited to 
DXLIM (fluid within boundary layer). 

ULlMI: Maximum permissible dimensionless 

velocity difference at inner free 
boundary (KIN = 2). 

ULIME: Maximum permissible dimensionless 

velocity difference at outer free 
boundary (KEX = 2). 



NF: 


Specifies the number of dependent 
variables to be solved for 



UE: 


TE: 


PE: 


Axial velocity of the outer fluid 
(m/s) 

Temperature of the outer fluid 
( K) 

o 

Pressure of the outer fluid (N/m ) 
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Data Cards. 


There are 120 data cards to be read in from subroutines 
THERM and OUTPUT. 

The first 99 cards include the values of the chemical 
equilibrium constants RC(60,NR) of the NR reactions and 
the enthalpies HT(60,NS) of the NS species. The data are 
read in at intervals of 100 K for temperatures up to 
6000 K, so for example, HT(35,2) refers to the enthalpy 
at 3500 K of species 2 (oxygen). 

The remaining 21 data cards contain alphanumeric data 
which are used to provide headings for the print-out. 
These are read in from OUTPUT and supply information on: 

i) the dependent variables of the calculation (10 cards) 

ii) the chemical reactions considered (NR cards) 

iii) the species present (NS cards) 
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